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ABSTRACT 

We present the results of a series of numerical simulations of compressible, self- 
gravitating hydrodynamic turbulence of cluster-forming clumps in molecular clouds. 
We examine the role that turbulence has in the formation of gravitationally bound 
cores, studying the dynamical state, internal structure and bulk properties of these 
cores. Complex structure in turbulent clumps is formed provided that the damping 
time of the turbulence, idamp is longer than the gravitational free-fall time ig in a re- 
gion. We find a variety of density and infall velocity structures among the cores in the 
simulation, including cores that resemble the Larson-Penston collapse of an isothermal 
sphere (p oc r~^) as well as cores that resemble the McLaughlin-Pudritz collapse of 
logatropic spheres {p ex r~^). The specific angular momentum profiles range between 
j oc r^ — r^. The masses of the bound cores that form are well-fit by the turbulent 
mass spectrum of iPadoan fc Nordlundl l)2002|) . while the specific angular momentum 
distribution can be fit by a broken power law. While our hydrodynamic simulations 
reproduce many of the observed properties of cores, we find an upper limit for the star 
formation efficiency (SFE) in clusters of 40-50 per cent. 
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1 INTRODUCTION 

Two of the most fundamental physical properties of stars 
that must be explained by a complete theory star forma- 
tion are the distribution of their initial masses (the initial 
mass function or IMF) and their initial angular momenta 
(the initial angular momentum function - or lAMF). The 
past decade has witnessed significant progress in measure- 
ment of the IMFs of stars in young star clusters. Infrared 
camera surveys of young stars within forming star clusters 
in molecular clouds reveal that the stellar mass distribution 
of stars is universal and foll ows the IMF of field stars (see, 
for example, the reviews of iMvers. Evans fc Ohashil 200C ; 



or binary stars will form) that is essentially the same as 
IMF (e.g.lMotte Andre fc Nerilll998l : iTesti fc Sargen3ll998l : 
I Johnstone et alJl200Cl) . This suggests the powerful and sim- 
ple hypothesis that the IMF is inherited from the mass spec- 
trum of molecular cloud cores. In such a picture, one-to- 
one correspondence between the IMF and the CMF natu- 
rally arises if the gravitational collapse of molecular cores 
into protostellar discs is followed by the efficient accretion 
of this material onto their central young stellar objects 
(YSOs). Such an accretio n scenario could in principl e ac- 
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Lada fc Lada 2003 



Molecular cloud cores in which individual (or binary) 
stars are born have been intensively studied for more than a 
decade, and complete surve ys of their physical properties 
have become available (e.g. Ijiiina. Myers fc Adama Il999l: 
[johnstono ct al. 2001). Millimetre and submillimetre sur- 
veys of the dense gas cores within cluster forming clumps 
have a core mass function or CMF (in which individual 
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count for both low mass JShu. Ad ams fc Lizano||l987 ) and 
high mass star formation (e.g. lMcLauehlin fc Pudritd 19971 : 
iMcKee fc Tan 2003) 

Excellent quantitative information is also available 
about the lAMF of young stellar objects. A recent study of 
254 stars in the Orion nebula as an example, finds that their 
YSO rotation periods perio ds are statistically consistent 
with a uniform distribution JStassun et alJll999^ ■ Molecu- 
lar cloud cores are observed to have specific angular mo- 
menta that are broadly distributed around a characteristic 
value of lO'^^ cm^ s~^ - much larger than the average angu- 
lar momenta of stars. This is the observational basis for the 
well-known "angular momentum problem" for young stel- 
lar objects - that the total amount of angular momentum 
seen in molecular cloud cores exceeds that measured for the 
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lAMF by several orders of magnitude. Since YSOs acquire 
their initial spin from their parental cores, it is important 
to understand how the core angular momentum function 
(CAMF) arises. 

These basic questions emphasize the need for theoretical 
and numerical work on the origin and evolution of molecu- 
lar cloud cores. A promising current idea is that molecu- 
lar cloud cores result from the turbulent fragmentation of 
clouds and their clumps. Cores are often observed to be ar- 
ranged in long, filamentary structures, a feature that can 
be qualitatively reproduced i n supersonic turbulence simula- 
tions a s a network of shocks ([ Porter^oxMUB^^^Woo^warc 



1994| : jBurkert fc Bodenheimeil_ l2000l: iKlessen fc Burkeri 



2000'; 'Klesscn k Burkcrt"2001'; "Ostr iker. Stone fc Gammie 
2001 ; Padoan ct al. 2001) . The statistical properties of the 



turbulence may naturally account for the distribution of 
masses and other properties of mo lecular cloud cores, such 
as their specific angular mome nta JBurkert fc Bodenheimeii 
l2000l:lKlessen fc Burkerl200lD 

This paper examines the origin, internal structure, and 
evolution of gravitationally bound structures that are pro- 
duced in 3D hydrodynamic simulations of self-gravitating 
gas. Our simulations follow the hydrodynamic evolution of 
an initially uniform density fluid that is perturbed by a spec- 
trum of velocity fiuctuations that damp out with time. 

Support against self-gravity in real molecular clouds in- 
volves a combination of turbulent and magnetic effects, as 
thermal and rotational support is negligible. In purely hy- 
drodynamical models for clouds, only the self-gravity and 
turbulence matters. Unless feedba ck processes from star for- 
mation (eg. bipolar outfiows - see iMatzner fc McKeelll999l) 
can replenish the turbulent energy within a dense clump in a 
molecular cloud, turbulent energy will damp out with time. 

We establish that the ratio tdamp/ift determines 
the time it takes until gravitational collapse occurs; if 
idamp/iff < 1 then turbulent fluctuations must damp before 
significant gravitational evolution can occur, and thus col- 
lapse is significantly delayed. If fdamp/iff > 1, the fiuid will 
experience gravitational collapse considerably earlier, while 
significant turbulent motions remain, such that there will be 
multiple gravitationally-bound cores. 

We find that turbulent fragmentation of a system of 
shocks not only accounts for the observed CMF - as other 
numerical simulations have shown - but for their angular 
momentum distribution (CAMF) as well, as the shocks will 
typically collide at an oblique angle. We also find that as 
long as clumps and cores remain turbulent, then their radial 
density profiles more closely resembl e those of turbulent lo- 
gatropes, analyzed bv McLaughlin fc Pudrita il996l) . than 
isothermal spheres. 

We discuss the technical aspects of our simulations in 
Section |5| where we also describe a new "watershed" core 
finding algorithm. In Section |3 we follow the energetics of 
our turbulent system, examining in particular the role of 
damping of turbulent energy and the applicability of the 
virial theorem and Bonnor-Ebert equilibria. We then go on 
to show the internal structures of cores (Section ^ as well 
as statistical properties of the ensembles of cores such as the 
CMF and CAMF (Section 01 . 



Run 


nj 


Mrms 


mTOT(ni0)* 


L(pc)* 


mj (niQ) 


A2 


1.1 


2.0 


12.5 


0.10 


11.6 


A5 


1.1 


5.0 


12.5 


0.10 


11.6 


B2 


4.6 


2.0 


105.1 


0.32 


22.9 


B5 


4.6 


5.0 


105.1 


0.32 


22.9 



Table 1 . The parameters for each of the runs that are highlighted 
in this paper. Each simulation was performed on 4 processors at 
a resolution of 256'^ grid points, with a kinetic energy spectrum 
described by a k~^^'^ power law for wavelengths shorter than 
one-quarter of the length of one side of the simulation box. 

The original clump parameters drawn from lLada. Bally &: Starld 
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Figure 1. Evolution of the total energy in our four simulations, 
normalized to the initial total energy of each run. The abscissa 
measures the time evolution in terms of the flow crossing time of 
each run. 



2 SIMULATIONS 

Our simulations are performed using t he parallel imple- 
menta tion of the ZEUS code (zEUS-MP) of lStone fc NormanI 
(Il992l). We us e the parallelized FFTW libraries 
JFrigQ fc Johnsonl |l993) to calculate the gravitational 
potential from the Poisson equation. The simulations 
were performed using a Compaq AlphaServer SC40 at 
the SHARCNET supercomputing facility at McMaster 
University. Our runs were performed at 256'^ resolution on 
4 processors. 

An initial fiuctuating velocity field was generated by 
constructing a scalar field with a power spectrum with of the 
form fc~"", where k = 27r/A is the wavenumber for a par- 
ticular isotropic Fourier mode corresponding to wavelength 
A, and n is the three-dimensional power-law index for the 
fluctuation spectrum. We use n = 11/3, the power spectrum 
representing Kolmogorov turbulence. For true turbulent mo- 
tion, the phases of these plane waves are correlated, but the 
details of this correlation are not well understood. We ap- 
ply initial random phases to these plane waves, and give 
each plane wave a randomly-oriented direction in the com- 
putational grid. We make the initial fluctuation spectrum 
divergence-free (V ■ v = 0) by taking the curl (or, in the 
momentum-space representation, taking the cross-product 
with the wavevector — ik) of this vector field. The net result 
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for the momentum-space velocity field is 



5vi 



= Ik 



,i/2 i4> 



e"*(-jkx u)/|kx u|' 



(1) 



where u is the random unit vector giving the direction of 
the plane wave. We also truncate our velocity spectrum at 
a maximum wavelength of Amax = L/4:, where L is the to- 
tal length of one side of our periodic box. This procedure 
strongly curtails any nonlinear effects that may arise from 
the periodic boundary conditions, such as large scale shock 
waves. Our initial velocity field is then transformed into 
position-space, where it is normalized to the RMS Mach 
number for that simulation. As we do not include any forc- 
ing of the velocity field during the simulation, the initial 
turbulent velocity field gradually decays away with time. 

Our initial density field was uniform, and normalized 
such that a specified number of Jeans masses nj = mtot/mj 
were on the box. The thermal pressure was determined via 
an isothermal equation of state. 

The ZEUS code evolves the hydrodynamic equations 



dp 



+ V ■ (pv) 



dt 
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(2) 

(3) 

(4) 



(The ZEUS code can use the full magnetohydrodynamic equa- 
tions, but for this paper we are not using magnetic fields, as 
we want to isolate the purely hydrodynamic effects. We will 
consider magnetic fields in a future paper.) 



2.1 Observational basis for initial conditions 

Our simulations explored a region of parameter space in 
(nj, A^RMs)- We performed many simulations but in this pa- 
per discuss a total of 4 different but representative simula- 
tions within this parameter space, listed in Table The 
values f or nj were chosen fr o m th e masses and sizes of 
cores in iLada. Ballv fc Stark! il991f) . These are cores de- 
tected in CS emission, corresponding to mean densities of 
10* — 10^ cm~^. Given their masses of 10 — 100 m© and 
radii of 0.1 — 0.3 pc, and assuming a thermal temperatures 
of 20 K (typical for molecular cloud cores), these clumps 
cont ain a few (1 — 10) initial thermal Jeans masses each. 
The iLada. Ballv fc Stark ( 1991) data set also contains typi- 
cal linewidths for these cores of 1 — 2 kms~^. For an isother- 
mal fluid at 20 K, for which the sound speed is ~ 0.4 kms~^, 
these linewidths correspond to flows of a few times the sound 
speed. We adopt the values of Mrms = i'rms /cs = (2, 5) for 
the RMS turbulent Mach number for the fluid. We focus our 
analysis on four particular combinations of nj and A/rms, 
namely (nj,A^RMs) = (1.1,2.0) (Run A2), (1.1,5.0) (Run 
A5), (4.6,2.0) (Run B2), and (4.6,5.0) (Run B5). 

While the calculations were performed in dimensionless 
units, they can be scaled to real values by choosing two of 
the total mass mtot , the total length of one side of the box 
L, or the temperature T, via 



me 



pc) V20Ky 



(5) 



It should be noted that Equation @ cannot realistically be 
scaled to arbitrarily large or small values of mtot or L, as 



real molecular clouds will only have nj in the range of 1 — 5 
for a certain range of L. 

Our simulations were allowed to run until the local 
Jeans length at some point within the computational grid 
was less than 4 pixels in length, a co ndition established t o 
avoid artificial fragmentation effects JTruelove et al.l[l998fl . 
For our simulations this works out to a density threshold 
that we can express in terms of the number of Jeans masses 
on the grid, and the initial or mean density: 

1/3 



Po 



80 



-2/3 



: 4148.6n 



-2/3 



(6) 



where A'^ is the number of pixels along one edge of the box, 
and po is the initial density of the fluid. 



2.2 Watershed algorithm for finding fluctuations 

We use a watershed algori thm (e.g. IVincent fc Soilldll99ll : 
iMangan fc Whitakerlll999r) to find fluctuations in our sim- 
ulations. We chose this algorithm because of its intuitive 
simplicity and speed of execution. As we show in Appendix 
1X1 it is c omparable to the commonly use d CLUMPFIND al- 
gorithm JWifliams. de Geus fc Blitzl^l994^ ■ in the limit of 
infinitesimal contour levels. The algorithm consists of the 
following procedure, iterated over the entire grid. 

At each initial cell, we find the local gradient vector 
and move to the next cell in that direction; we repeat this 
process until we reach a local maximum. At each cell that 
we pass through in this path, we assign it an identification 
number q for this path. If we come to a cell that already has 
an identification number q' , we re-assign all of the cells on 
the current path q to the identification number of the new 
path q' . This way, all cells that are associated with a single 
local maximum are assigned to the same fluctuation. 

Our algorithm is extremely sensitive to small fluctua- 
tions in density, resulting in a large number of fluctuations 
found at the numerical resolution limit. We therefore flrst 
smooth the data by averaging each cell with its nearest 
neighbours in each direction before running the fluctuation- 
finding algorithm. This greatly reduces the number of fluc- 
tuations found with a radius of only a few pixels. After the 
fluctuation-finding routine had been run, we used the origi- 
nal, unsmoothed data to calculate the fluctuation properties. 

When all of the grid cells have been assigned to fluc- 
tuations, we can look at the boundaries of each fluctuation 
(defined as a cell which has a 6-neighbour that belongs to 
a different fiuctuation), and calculate the surface density 
(from which we can extract the surface pressure), and the 
mean radius of the fluctuation. We can also calculate the 
mass, mean velocity and internal velocity dispersion of the 
fluctuations from the fluctuation data. 

We show in Appendix ^ that our watershed algo- 
rithm gives results that agree with the results returned 
by CLUMPFIND when applied to the same numerical data. 
Both the CLUMPFIND algorithm and our watershed algo- 
rithm flnd a similar number of bound, self-gravitating fluc- 
tuations that we define as "cores" (see Section r^.H ; however, 
the watershed algorithm finds many more unbound density 
fluctuations than CLUMPFIND. As the density resolution of 
CLUMPFIND is increased, the number of unbound fluctuations 
increases, approaching the number found by the watershed 
algorithm. 
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10 



RMS Mach Number 



Figure 2. The ratio of tdamp/iff predicted by Equation J7J- The 
solid lines, from right to left, represent ratios of 1.5,2.0,3.0 and 
4.0. 



As one of the goals of this paper is to examine what 
properties are needed for a density fragment to undergo 
gravitational collapse, we want to have a complete census of 
turbulent fluctuations. The watershed algorithm presented 
here is more sensitive to these unbound fluctuations, so we 
consider this definition of a fluctuation to be better suited 
for our purposes than the definition provided by CLUMPFIND. 



3 ENERGETICS 

Core formation within our simulations is driven by the ener- 
getics of self-gravity and turbulent decay. We plot the evo- 
lution of total energy Et = Skin + -Egrav in Fig. As our 
simulations employ an isothermal equation of state, total 
energy is not conserved (in a real system, this would corre- 
spond to energy losses via radiative cooling). The evolution 
of the systems can be divided into three phases: an initial 
phase that is dominated by shocks forming in reaction to 
the initial conditions; an intermediate phase that primarily 
involves the dissipation of kinetic energy in the shocks; and 
a late phase that is dominated by gravitational contraction 
of cores. 

The flrst phase lasts ~ 0.1 tfiow in each of our runs, 
marked in Fig. by the rollover at early times. 

Turbulent kinetic energy decays via the 
power-law i?K o c t ~^ i n a compressible fluid 
jMagJjO^ ^9gj |Ostriker^Gammte&StoM 



19991: IChristensson. Hi ndmarsh fc Brandenburd 1200 L 
Ostriker. Stone fc Gamniial200l|) . In a linear-log plot like 
Fig. this power-law relationship appears as a decaying 
exponential. A signiflcant portion of runs B5, A2, and 
especially A5 is spent in this phase. In run B2, the ki- 
netic energy damped sufficiently during the shock-forming 
phase, allowing gravitational collapse to proceed without a 
signiflcant decay phase. 

The final phase of these simulations begins when the 
turbulent energy has damped sufflciently that gravitational 
contraction can overcome turbulent pressure to force at least 
one core to collapse. In Fig.0 this is marked by a sharp de- 
crease in the total energy as the gravitational energy grows 



to large negative values. The Mach 5 runs have a larger ki- 
netic energy for the same amount of mass as the Mach 2 
runs, and so can resist gravitational collapse for more flow 
crossing times. Similarly, the B runs have more mass than 
the A runs, and so less decay of turbulent energy is needed 
before gravitational effects become dominant. 

The simulations with 1.1 Jeans masses take many flow 
crossing times before the turbulent motions damp suffi- 
ciently that self-gravity can overcome the combination of 
turbulent and thermal internal pressure; with turbulence 
having damped away, a significant amount of the fluid is 
gathered into a few large cores. In the simulation with 4.6 
Jeans masses, gravitational collapse occurs after only a few 
flow crossing times, so that the turbulence has not damped 
as much. While the flrst object to experience gravitational 
runaway collapse causes the simulation timesteps to become 
very small, thus preventing us from observing the evolution 
of any other objects, we still flnd many gravitationally bound 
cores amongst the wispy fllamentary structures in the fluid; 
we would thus expect to flnd a number of star-forming cores 
in such a fluid. 

Insight into the relationship between the time it takes 
for runaway gravitational collapse to occur can be found in 
the relationship between the turbulent damping time and 
the local gravitational free-fall time of an individual frag- 
ment. The damping time of the turbulence in the simulation 
is related to the flow crossing time, idamp = L/u = L/{csM), 
where L is the length of one side of the simulation box, 
u — CsM is the root-mean-square velocity of waves in the 
fluid, and M is the Mach number. As the gravitational col- 
lapse occurs in the enhanced density behind the shocks, 
the local free-fall time can be calculated from the post- 
shock density, tg = (32G'pi/37r)~^". For a system of strong 
isothermal shocks, pi ~ poM, where po is the average, pre- 
shock density of the fluid. Using Equation ^, the free-fall 
time can be expressed in terms of our nj parameter, and we 
can calculate the ratio of the damping time to the free-fall 
time to be 



ida 



ts 



= 3.25 nJ/^M" 



-1/2 



(7) 



If idamp/iff <C 1, we expect turbulence to damp away be- 
fore collapse can occur to any signiflcant degree, leading 
to a few main cores that dominate. For tdamp/iff 3> 1, 
we expect the fluid to collapse while turbulent motions re- 
main strong, thus resulting in many small cores located 
within the shocks. This limit cor responds to the standard 
turbulent fragmentation scenario JKlessen fc B urkerlil200CI : 
lOstriker. Stone fc GammidbOOlt IPadoan e t alj|200j) . 

We plot the relationship in Equation JTJ , along with the 
values for our simulations of (nj, Mrms), in Fig.|5| All of our 
runs have 1.5 < tdamp/iff < 4.0. The ratios for tdamp/iff for 
the runs are consistent with the time it takes for collapse to 
occur in Fig. in that the simulations with a larger ratio 
of idamp/tff collapse signiflcantly quicker. 



3.1 Virial equation 

We can gain a more detailed insight into the dynamics of 
these fluctuatio ns by analyzing the ter ms in the Eulerian 
virial equation jMcKee fc Zweibellll993) : 
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Figure 3. The importance of the different terms in the virial equation. S is the surface terms, W is the gravitational term, K is the 
kinetic energy term, and U is the internal energy term. The solid line corresponds to virial equilibrium, I' = 0. In this figure and the 
following figures, run A2 is plotted in the upper-left, A5 in the upper-right, B2 in the lower-left, and B5 in the lower-right. 



1 •• Id 

2 ^ 2di 

where 



/ = 



{pvr 



pr dV 



dS = U + K + W + S = -i' (8) 



(9) 



(7 = 3/ PdV 

Iv 

K = I pv^dV 
'v 



W 



pr ■ V$dV 



[Pr + r ■ (pw)] ■ dS 



(10) 

(11) 
(12) 
(13) 



where U is the thermal energy of the fluctuation, K is the 
internal kinetic energy (including both collapse motion and 
internal turbulence), and W is the gravitational potential 
energy. Each dot over a symbol represents a time derivative. 
The terms U and K act to support the fluctuation, while 
W drives gravitational collapse. The two terms in S rep- 
resent the thermal pressure and turbulent pressure on the 
surface of the cloud. The second term on the left-hand side 
of Equation ||SJ| is a time derivative of the moment-of-inertia 



flux through the surface of the fluctuation. We are unable 
to directly calculate this quantity, even though it could be 
non- negligible and of either sign. We will treat it as if it were 
negligible, such that we will call 7' = virial equilibrium. It 
should be remembered that this will be merely an estimate 
of the true virial equilibrium of the fluctuations. 

The internal energy and kinetic energy terms must al- 
ways be positive. For an isolated fluctuation, the gravita- 
tional term will be negative (counteracting the effects of 
thermal pressure), but the fluctuations in our simulations 
are not in isolation — they can be disrupted by tidal eflects 
from nearby fluctuations. On average, about a quarter of the 
fluctuations found in each of our simulations have a gravi- 
tational term which acts to increase I. The surface pressure 
term will usually be negative (but in theory could be pos- 
itive); for all of the fluctuations in each of our simulations, 
the surface pressure acts as a confining force. 

We compare the surface terms in the virial equation to 
the gravitational term in Fig.|3] For both axes, we normalize 
the virial terms to the internal kinetic and thermal terms of 
the virial equation. The solid line marks the locus of virial 
equilibrium on the graph. The fluctuations that are the most 
out of equilibrium (as measure by the largest values of | J'|) 
cluster at the most negative values of W/{K + U) (as a 



6 Tilley, D. A. and Pudritz, R. E. 




10- 



10"'' 10-2 10° 

m/mBE-Crit 



102 




TQ-6 1Q-4 1Q-2 100 102 

m/mBE-Crit 




10-= lO-'^ 10-^ 10-2 10-"" 10° 10^ 

m/mBE-Crit 



10^ 








1 


1 


■ 


10° 


- 










- 


0-2 


- 








P 


- 


0-4 


i 


^ 


(« 






- 


0-6 




, 


o 


o 






10 


-6 


10 


-4 


m/mBE- 


10° 

Crit 


U 



Figure 4. Comparison of tlie stability from the virial equation, versus the stability from the Bonnor-Ebert criterion. The diamonds 
represent fluctuations with /' > 0; the stars represent cores with I' < 0. 



strongly self-gravitating core is expected to be) and the least 
negative values of S/{K + U). 

Most of the fluctuations in each of the runs have 7' > 
(for the plots in Fig. |3 this is equivalent to the fluctuation 
lying above the line I' — 0), indicating that they are likely 
to disperse. These fluctuations have relatively little gravita- 
tional energy, but the surface pressure can be signiflcant. In 
the A runs, only a tiny number of fluctuations are out of 
virial equilibrium with /' < 0, such that they will collapse. 
Even these collapse-unstable fluctuations are not far out of 
equilibrium, lying close to the virial equilibrium line. We see 
a different behaviour in the B runs; there are many fluctu- 
ations that are virial-coUapse-unstable, with some fluctua- 
tions that are far out of equilibrium. This is not surprising, 
as the A runs have had time for turbulence to damp away 
and for the system to reach a quasi-equilibrium state, but 
we were forced by our Jeans resolution criterion to stop the 
evolution of the B runs while significant turbulent motions 



We use the virial equation to divide our set of density 
fluctuations into two subsets. The first is what we will call 
bound cores . These cores have 7 < 0, indicating that they 
are bound by pressure and/or gravity. Every density fluctu- 
ation that is not a bound core will be called an unbound 
fluctuation. 



3.2 Bonnor-Ebert stability 

For each fluctuation we can calculate the mass of a criti- 
cal Bonnor-Ebert sphere from the velocity dispersion and 
surface pressure |Ebert 195f^) : 



m-BE 



17.6 
AttG^Ps 



(14) 



where Ps is the surface pressure, Ca is the sound speed, and 
G is the gravitational constant, tube represents the maxi- 
mum mass the fluctuation could have, given its temperature 
and surface pressure, before it becomes unstable to collapse. 
If the mass of an individual fluctuation is greater than tube, 
it suggests that that fluctuation might be in a state of col- 
lapse. Conversely, if the mass of a fluctuation is less than 
rriBE, the fluctuation is either stable or not gravitationally 
bound. 

Out of the few thousand objects our algorithm identi- 
fies as fiuctuations, only the most massive one or two fluctu- 
ations are Bonnor-Ebert-unstable. The surface pressures of 
the other fluctuations are low enough that the Bonnor-Ebert 
critical mass is larger than the actual fluctuation mass. 

We can compare the stability criterion established by 
Bonnor and Ebert to the stability criterion from the virial 
equation. We plot the absolute value of 7' against the ratio 
of the mass of the fluctuation to the Bonnor-Ebert critical 
mass in Fig.0] the fluctuations marked by stars are the ones 
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Figure 5. The ratio of central density to surface density, as a function of m/r?iBE- The horizontal line in each plot is Pc/ps = 2. 



that are likely to collapse from virial arguments. We can 
see that the fluctuations that are Bonnor-Ebert-unstable are 
also the most virial-unstable, although many fluctuations 
which are stable according to the Bonnor-Ebert criterion 
are unstable according to the virial equation. We also note 
that the fluctuations with /' < (marked by stars) follow a 
distribution with a steeper slope than the fluctuations with 
/' > (marked by diamonds). 

Our fluctuation-flnding method is sensitive to any local 
density fluctuation, but not all density maxima in our sim- 
ulations would be detected in surveys for cores. In observed 
maps of molecular clouds, cores are detected by virtue of 
having a substantial density contrast between centre and 
edge. We examine the relationship between this density 
contrast and m/meE in Fig. |^ A substantial fraction of 
fluctuations in the B runs have a significant density con- 
trast between surface and center, but only a handful of 
fluctuations in the A runs have Pc/ps > 2. The region 
Pel Ps > 2 is the t hree-dimensional analogue of what the 
Ijiiina. Mvers fc AdamSi (.1999.') data set would consider a 
core. 



3.3 Virial masses 

It is useful to compare the true mass of the fluctuations 
found in our simulations to the mass that would be predicted 
by standard observation techniques. A common method used 



with molecular line observations is to assume a core is in 
virial equilibrium (such that 2K -I- W = 0; this doesn't con- 
sider surface effects), and calculate the appropriate mass 
given the size and velocity dispersion of the core. If one as- 
sumes that the cores are un i form s pheres (as if often done, 
for example iGoodman et alj il993l l: ljiiina. Mvers fc Adamj 
L1999:) 1. then the virial mass is rrivir = f ^tj^i where a is the 
turbulent velocity dispersion. We plot the relationship be- 
tween this predicted virial mass and the true mass of our 
fluctuations in Fig. |S| For all four runs, this measure of 
the virial mass generally overestimates the mass of the fluc- 
tuations by two orders of magnitude; the average ratio of 
m^ir/m is 40.6 for run A2, 66.4 for run A5, 211.9 for run 
B2, and 178.2 for run B5. Furthermore, we find that the re- 
lationship between nivir and m is approximately a straight 
line on a log-log plot, but with a slope less than one (0.53 
for runs A2 and A5, 0.37 for run B2, and 0.45 for run B5). 

Observational comparisons between the mass estimated 
by virial techniques and masses estimated from other meth- 
ods also find that the assumption of virial equilibrium leads 
to masses l arger than that m e asured by other methods. 
A study by Ivan der Tak et alj ll2000h found that rrivir ~ 
2.77mother, and that the slope of the log(mvir) — log(mothcr) 
plot to be slightly less than 1.0. Our results suggest a much 
larger discrepancy between true mass and virial-estimated 
mass. 

The discrepancy between TTivir and m is likely due to 
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Figure 6. The relationship between the mass that would be calculated by applying the virial theorem and the actual mass of the 
fluctuation. The solid line is the equality rrXvir = va. 



three effects. First, many of the fluctuations in our simula- 
tions are not in equilibrium; turbulent motions play a large 
role in the dynamics of these fluctuations. This is especially 
true in the B runs, which have more Jeans masses on the 
computational grid. Second, the simple virial mass estima- 
tor neglects the effects of the surface terms, which as we 
have already shown are dynamically signiflcant, especially 
for low-mass fluctuations. Finally, for cores that are cen- 
trally condensed, the gravitational energy will be greater 
than that of a uniform core. This will reduce the virial mass 
estimate, as less mass is needed for there to be virial equi- 
librium between kinetic and gravitational energies. 



4 PHYSICAL PROPERTIES OF SIMULATED 
CORES 

4.1 Global images and line maps 

A volume rendering image of a slab of our data cube from 
one of the runs we discuss (B5) is presented in Fig. |7| This 
particular image occurs at 0.5 flow-crossing times, and is not 
centered on any particular feature as Figs. ISl— 1101 are. As 
our fluctuation-finding algorithm associates any local den- 
sity peak with an individual fiuctuation, it is not surprising 
that we find many fluctuations in our simulations. The fluc- 



tuations are typically arranged along the wispy fllaments 
that are the result of interacting shocks. 

Observations of molecular cloud cores are unable to re- 
solve the three-dimensional structures of the cores, as we 
are able to here. Instead, they must resort to spectral line 
information, inferring the structure from the shape of the 
spectral lines. We provide for comparison in Fig.|H|line maps 
of our simulations at their final time outputs, overlayed on 
simulated observations of the column density. There are sig- 
nificant differences between the evolution of the A runs with 
1.1 initial Jeans masses on the grid and the B runs with 4.6 
initial Jeans masses. There does not appear to be a signif- 
icant difference between runs with different initial kinetic 
energies that have the same number of Jeans masses. 

The line maps are generated by treating the fiuid as if it 
were in the optically thin limit. For densities of ~ 10^ cm-"^ 
as found in our simulations, this might correspond to a tracer 
like CS. We take lines-of-sight along one of the axes of the 
grid, and for each lin e-of-sight, we calculate the line pr ofile 
via the equation (e.g. IWalker. Naravanan fc Bosslll994) 



I{Vohs 



Epe 



-(«o 



-^iu.YI< 



(15) 



where Uobs is the velocity channel being observed, wios is the 
true line-of-sight velocity in that cell, and the sum is over 
all of the pixels along the line-of-sight. Equation 11511 as- 
sumes that the line is optically thin, so that the details of 
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Figure 7. Volume rendering of the density in a 256x256x32 slab from one of our simulation runs, B5. 
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Figure 8. Line profiles and contour maps of tiie column density of our simulations. The abscissa on each of the hne profiles is in units 
of the Mach number. The line profiles in each 32x32 pixel block are binned together; the column density map has been convolved with 
a Gaussian of FWHM of 32 pixels as well, to simulate a typical observation. Note that the data has been shifted so that the peak in the 
column density is located at the centre of the image. The contours are at (10,20,30,40,50,60,70,80,90) per cent of the maximum column 
density value. 



line emission depend only on the amount of matter present 
to radiate. The Gaussian term includes the effects of ther- 
mal broadening. The resulting set of line profiles are then 
normalized with respect to the maximum line intensity in 
the map. The lines-of-sight are then averaged into bins in 
X and y, and displayed over the column density map. The 
area covered by the line profiles on the column density map 
indicates the size and location of each bin (generally 32x32 
pixels) . 

The contour map in Fig. |5] is generated by convolving 
the integrated column density with a Gaussian filter, with a 
FWHM of 32 pixels. The contour levels are linearly spaced in 
column density, with values of (10,20,30,40,50,60,70,80,90) 
per cent of the maximum density value. The combination 
of the convolved column density map and line map simu- 
lates a typical observation of a molecular cloud clump, where 
the resolution of the data has been reduced by the resolu- 
tion of the telescope, and the structure smeared out by the 
beam (for comparison, the unconvolved column density map 
can be found in Fig. |^ with logarithmic contours). To put 
these numbers on an observational basis, at a distance of 



1000 pc, this "beamwidth" would correspond to an angular 
beamwidth of 2.6 arcsec for the A runs, and 8.3 arcsec for 
the B runs, if we use the size scaling found in Table (for 
a cloud like Orion at ~ 500 pc, these numbers work out to 

5.2 arcsec for the A runs, and 16.6 arcsec for the B runs). 
The Very Large Array, for comparison, has a FWHM of 8.4 
arcsec at 3.6 cm wavelengths in D-array, and 2.8 arcsec at 

1.3 cm wavelengths. Thus, the maps in Fig.|5]are analogous 
to the type of line map produced by these radio telescopes. 

In the lower- Jeans-mass runs, the initial shocks created 
by the turbulence can not accumulate enough mass for grav- 
ity to overcome thermal and turbulent pressure. As a con- 
sequence, the shocks must dissipate before gravity becomes 
significant. Gravitational collapse only occurs once the tur- 
bulent energy can decay sufficiently for the self-gravity to 
overcome the turbulent pressure. Only a few core-like struc- 
tures can be identified, as the fiuid motions are dominated 
by the self-gravity of the main core, rather than being bro- 
ken up by turbulent motions. The line profiles in this map 
tend to be smooth and symmetric, although the line profile 
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Figure 9. Unconvolved contour maps of the column density of our simulations. As in Fig.|S] we have shifted the data so that the peak of 
the column density distribution is in the centre of the image. The contour levels are logarithmically spaced, starting at a column density 
of 2.4 X 10^^ cm~'^ for the A runs (when scaled to the 0.1 pc size in Table HI and a column density of 1.9 X 10^^ cm~^ for the B runs 
(when scaled to the 0.32 pc size in Tabic IH . and increasing a factor of 1.5 with each successive contour level. 



through the centre of the large core for both of the A runs 
is quite wide. 

In the higher- Jeans-mass simulations, as in the low- 
mass runs, the initial shocks are not strong enough to in- 
duce gravitational collapse. However, it takes significantly 
less time for the turbulent energies to decay to the point 
where gravitational collapse can occur, because tg is shorter. 
As a result, significant turbulent motions still remain, and 
the fluid is broken up into many fragments. We can thus 
identify many potential cores in these simulations. The line 
proflles are significantly more asymmetric for the B runs 
than the A runs. This refiects the importance of turbulent 
motions on the line profiles. 

Due to the stochastic nature of the turbulence, some 
dense core will collapse before any other fluctuation can 
evolve to a signiflcant degree. As our criterion for stopping 
the simulations is based on a density threshold, the simu- 
lation will end before many regions of the simulation can 
collapse. This is not a problem in runs A, as the turbu- 
lent effects have damped out and a single core dominates 
the simulation. Stochastic effects are a concern in runs B, as 
turbulent fluctuations may cause additional cores to undergo 



runaway growth in the same manner as the first significant 
core that forms. This could be a potential bias when inter- 
preting statistical information about the cores, such as the 
mass distributions. 



4.2 Internal structure of cores 

The fiuctuations and cores in our simulations can have com- 
plicated structures. We show slices through one plane of the 
largest core in each of the runs in Fig. llUI that show both the 
density structure and velocity field in that slice. The A runs 
exhibit the fairly symmetric structure of the largest core 
that dominates the image. The B runs, however, show an 
intricate density and velocity field as the fiuid is channeled 
into the larger cores primarily along filamentary structures. 
In the absence of turbulence, a pressure-bounded 
isothermal sphere that is gravitationally unstable will col- 
lapse. The details are dev eloped analytically by iLarsonl 
(1969) and IPenstorJ il969h . and have a self-similar form. 
The Larson/Penston solution consists of a flat inner region, 
transitioning to a r~^ profile on the outside. The collapse 
proceeds asymptotically towards a density singularity (in a 
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Figure 10. A slice tlirougli tlie centre of tlie most massive core in each of tiie simulation runs. The greyscale contours show the density 
distribution; the arrows represent the velocity flow in the plane of the image. The length of the velocity vectors are in proportion to the 
speed of the fluid flow, normalized such that a vector with a length equal to the distance between the tails of the vectors has a speed of 
Mach 2. 



real system, collapse would be halted by opacity and ther- 
monuclear fusion as a protostar). The infall velocity reaches 
a maximum of 3.3 times the sound speed just outside of 
the flat central density r egion. Numerical calculations by 
iFoster fc Chevalieii (|l993i) confirm this behaviour. 

The presence of turbulent motions complicates this 
analysis. A phenomenological attempt to explain turbu- 
lent linewidths is a logatropic equation of state, P/Pq = 
1 -I- Alnp/po, which has the effect of providing relatively 



more pressure support in low density regions than in 
the high-density regions, in comparison to an isothermal 
equation of state [A - - ^ 0.2 is an adjustable parameter). 
iMcLaughlin fc Pudrita (119971) performed the analogous cal- 
culation to the Larson-Penston self-similar solution, using 
this logatropic equation of state rather than the isother- 
mal equation of state. As with the Larson-Penston solu- 
tion, the logatropic collapse consists of a flat inner part, 
but the outer parts follow a r~^ profile. Three-dimensional 
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Figure 11. Density and infall velocity plots for run A2. The most massive core is on the top, followed by the second and third most 
massive. The first core is bound according to the criterion in Section |3. II the second and third cores are both unbound. The abscissa is 
in units of pixels; for a box size of 0.1 pc (from Table HI . one pixel corresponds to 81 AU. 



collapse calculations of a logatrgpic fl uid reproduce these 
results jReid. Pudritz fc Wadslgvli2002l) . 

We plot the density and velocity profiles of the three 
most massive cores in each simulation in Figs. I11I14I These 
plots are generated by finding the average value of the den- 
sity and radial velocity within each shell of radius r; the 
vertical bars indicate the total range of values that the den- 
sity or radial velocity has at that particular radius. The most 
massive cores in each run have a p oc r~^ profile everywhere 
(corresponding to the singularity formation of the Larson- 



Penston collapse). The lower- mass cores can have consider- 
ably shallower density profiles that could indicate that the 
core is at an earlier stage in the Larson-Penston solution, 
when only the flat central part of the self-similar solution is 
contained within the truncation radius, or that the core is 
behaving as the logatropic McLaughlin-Pudritz solution. 

The degree of asymmetry of each core at a given density 
level is reflected by the scatter at that density. We can see 
that this scatter is typically a factor of ~ 5 in radius, with 
most of the points within ~ 3. This implies that the cores in 
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Figure 12. Density and infall velocity plots for run A5. The most massive core is on the top, followed by the second and third most 
massive. The first and second cores are bound according to the criterion in Section lij.il the third core is unbound. The abscissa is in 
units of pixels; for a box size of 0.1 pc (from Table HI . one pixel corresponds to 81 AU. 



our simulation are elongated. The largest cores in our B runs 
are less spherical then the most massive cores in the A runs 
- as turbulence is significantly damped in the A runs at this 
time, this is a natural result of the balance of self-gravity 
and thermal pressure (as opposed to the B runs, which are 
still closely tied to the planar nature of the parent turbulent 
flow) 

Our results indicate that the more developed cores have 
structure that is dominated by gravitational effects - they 
exhibit the r~^ density profile of a collapsing isothermal 



sphere. The less- developed cores and unbounded fiuctua- 
tions have significant non-thermal support (remember, we 
are forced to stop when the most advanced core violates 
the lTruelove et alJ il997ft criterion), and so exhibit density 
profiles that can be much shallower, resembling that of a 
logatropic sphere or shallower. This reinforces the idea that 
these fluctuations will collapse only when their internal tur- 
bulent motions are sufficiently damped in order that gravity 
can dominate. 

The behaviour of the infall velocity profiles indicate 
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Figure 13. Density and infall velocity plots for run B2. The most massive core is on the top, followed by the second and third most 
massive; all three cores are bound according to the criterion in Section |3. II The abscissa is in units of pixels; for a box size of 0.32 pc 
(from Table m. one pixel corresponds to 258 AU. 



that, whatever their shape, the cores are not collapsing 
spherically - in some directions, they are not collapsing at 
all. This results in the large spread of Vr at any given radius. 
We will also note that a few of the lower-mass cores in each 
of our runs appear to have no infall at all in their centres. 
This suggests that their interiors are in hydrostatic balance, 
at least temporarily. 

We examine the rotational behaviour of fluctuations 
and cores by plotting the magnitude of the specific angular 
momentum j = J/M (where J is the total angular momen- 



tum of a parcel of fluid and M is that parcel's mass) as a 
function of radius in Figs. I15I16I As with the radial veloc- 
ities, the specific angular momentum profile consists of an 
envelope, with a range of values at smaller j. If the core were 
to undergo bulk rotation, a similar distribution of |j| would 
be expected, as the fluid elements in the plane of rotation 
will have maximal values of |j|, while the elements along 
the axis of rotation will have j ~ 0. The envelopes typi- 
cally follow & j (X r^ profile, with rj between 1 and 2. The 
most massive cores tend to have a central region of ~ 10 
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Figure 14. Density and infall velocity plots for run B5. The most massive core is on the top, followed by the second and third most 
massive; all three cores are bound according to the criterion in Section |3. II The abscissa is in units of pixels; for a box size of 0.32 pc 
(from Table m. one pixel corresponds to 258 AU. 



pixels with j ex r^ , with a slightly flatter profile in the 
outer regions of r^. This r^ profile is consistent with fluid 
elements in orbit around a p oc r~^ density profile (since 
Vrjy ~ constant), as we see in the more massive cores. As the 
smaller cores tend to have more turbulent support in their 
envelopes, the density profile is shallower, as shown previ- 
ously. As j ~ Vrotr ~ r^ p^'^ , this allows the fluid to support 
a speciflc angular momentum profile that can be as steep as 
r^'^ to r"^, as can be seen in Fig. 1151 

The maximum magnitude of j is ^ CsL/4, where Cs is 



the sound speed and L is the size of our simulation box (I//4 
is the maximum scale of the turbulence, which by virtue of 
our turbulent spectrum, is also the scale that contains the 
most turbulent kinetic energy). This is roughly the order- 
of-magnitude one would expect for cores collapsing out of 
oblique shocks in a turbulent fluid (the thickness of the shock 
r ~ L/AM, where M is the RMS Mach number of the fluid, 
is roughly the same size as the core that forms out of the 
shock; the velocity transverse to the shock is going to be on 
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Figure 15. Specific angular momentum profiles of the three most massive cores in each run. The cores plotted in the left-hand column 
are from Run A2, the cores plotted in the right-hand column are from Run A5. The abscissa is in units of pixels; for a box size of 0.1 pc 
(from Table0, one pixel corresponds to 81 AU. 



the order of ii ~ CsM, resulting in |j| r-^ vr '~^ CsL/4. Note 
that this is independent of M) . 



5 STATISTICAL PROPERTIES OF CORES 

5.1 Mass distribution 

The mass distribution for the fluctuations in our runs are 
plotted in Fig. 1171 We have plotted the entire data set of 



fluctuations for each of our run with a dotted line, and the 
subset of bound cores (see Section 13.11 with a solid line. 
We are using the differential mass distribution, such that 
dN (m) / dlog{m) counts the number of fluctuations with 
mass between log(7n) and log(m) + dlog(7n). 

The entire fluctuation set results in a distribution that 
is similar to a log-normal form at the peak mass and lower, 
but with a Salpeter power-law form at higher masses. The 
bound core subset consists of most of the massive cores, with 
a peak at '^ 2 x 10~^mj. As an example, for a clump with 
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Figure 16. Specific angular momentum profiles of the three most massive cores in each run. The cores plotted in the left-hand column 
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a size L = 0.32 pc, T = 20 K, nj = 4.6 (from Table [3, 
Equation @ states that witot = 105.1 tuq, so that the 
peak in the IMF is at ~ 0.05 niQ. Above this peak, the 
distribution is a power law with index ~ —1.3 for the B 
runs (the Salpeter value for the IMF of stars is —1.35); the 
A runs only have a few cores that are bound. While there 
are not sufficient bound cores in the A runs to perform a 
meaningful statistical analysis, we note that observations of 
isolated star- forming regions such as Taurus can have IMFs 
that are approximately flat (iLuhman et al.ll2003f) . 



The fits to the mass spectrum in Fig. 1171 ar e calculated 
using the analysis of lPadoan fc Nordlundl l|2003)- Their cal- 
culation posits that the cores out of which stars form are 
created as the result of the fragmentation of shocks. In this 
scenario, the high-mass slope a = —3/(4 — /3) is determined 
by turbulent power spectrum /3 of the fluid. Here, /3 is index 
of the one-dimensional power spectrum, calculated by inte- 
grating the three-dimensional power spectrum over a sphere 
in fc-space. The three-dimensional power spectra of our sim- 
ulations are given by {Sw^Y', where (5vfc is given by Equa- 
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tion Q. Thus, P — n — 2, where n is our three-dimensional 
power-law index from Equation Q. For n = 11/3, /3 = 5/3 
and a = —9/7 — —1.29. In the Padoan & Nordlund formal- 
ism, the width of the mass spectrum is determined primar- 
ily by the turbulent Mach number; the theoretical curves 
we present in Fig. 1171 are calculated using the initial RMS 
Mach number of each simulation. We scale this theoretical 
mass distribution such that the total number of cores are the 
same as the measured distribution for the bound cores from 
the simulation, and that the total mass contained within 
cores with mass less than 0.1m j are equal. The width of the 
theoretical distribution and the high-mass slope match the 
simulation mass distributions quite well. 



5.2 Star formation efflciency 

The fraction of the total mass contained within the bound 
subset of cores is quite large - 75 per cent for run A2, 67 
per cent for run A5, 49 per cent for run B2, and 42 per 
cent for run B5. In order to turn this into a star formation 
efficiency, we also need to consider that not all of the mass 
within these cores will necessarily end up within the star - 
it is not clear what fraction of the core mass will actually be 



accreted. This is something that our isothermal equation of 
state fails to take into account, as the gas in our simulations 
can cool quicker than an ideal gas can. If most of the mass 
in our bound cores goes into stars then we have found an 
upper limit on the SFE of a clump. 

A lower limit for the efficiency of turning gas into stars 
in molecular cloud clumps (the star formation efficiency, or 
SFE) is esti mated to be 10 — 3 per cent using infrared 
observations (ILada fc Ladal2003r . The fraction of gas within 
the bound cores of our B simulations is not significantly 
greater than this. As outflows and stellar winds remove some 
of the mass of these cores as a star forms within, the resulting 
SFE of our B simulations will agree with the observed values 
for clustered star-forming regions. The A runs would have a 
SFE significantly higher than this, although they might be 
appropriate for isolated star-forming regions. 

5.3 Specific angular momentum distribution 

In Fig. 1181 we plot the distribution of fluctuation mean spe- 
cific angular momenta in the same manner as we did for the 
mass distribution. dN{j)/d\og{j) is the number of fluctua- 
tions with specific angular momentum between log(j) and 
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log(j)+dlog(j). Here, the specific angular momentum is de- 
fined as j = J/m, where J is the total angular momentum 
of the fluctuation and m is the mass of the fluctuation. 

The peak of this distribution for the entire data set is at 
^ CsL/A (For sound speeds of 0.4 kms~^ and clump length 
scales of 0.1 pc, this works out to 3.1 x lO'^^ cm^s~^). As 
we have shown in Section r4.2l this is the order-of-magnitude 
that one would expect from a network of oblique shocks. 
For the bound cores, the peak occurs at ~ bCsL. The over- 
all distribution appears to be symmetric, with a width of 
approximately one order-of-magnitude. The bound core dis- 
tribution is more asymmetric, with slightly more cores at 
lower values of j. We can fit the bound core distributions 
with a broken power law for the B runs, which qualitatively 
gives a good fit at both low-j (with a slope of ~ 1) and 
high-j (with a slope of ~ —2) values. 



6 DISCUSSION AND CONCLUSIONS 

Our set of simulations support the idea that molecular cloud 
cores, in which stars can form, arise from the turbulent frag- 



mentation of the parent cloud. We posit that the cores are 
forming out of oblique shocks that impart angular momen- 
tum, and thus a net rotation, to these cores. We have demon- 
strated that the evolution of the cores that form from molec- 
ular cloud turbulence is driven by the degree of turbulence 
within the cores. The mass distribution of cores that we find 
ap pears to be well-fit by the turbulent mass distribution 
of IPadoan fc Nordlundl ll2002r . suggesting that these cores 
were initially formed from shocks. Gravitational collapse can 
only occur when the turbulent energy has damped suffi- 
ciently. There are also many cores that are bound through 
turbulent surface pressure, although these are not in a state 
of dynamical collapse. The cores in general are not in virial 
equilibrium; indeed, many are quite far from virial equilib- 
rium. The most massive cores are supercritical in the sense 
of Bonnor and Ebert. 

Our simulations show that structure formation in a tur- 
bulent, self-gravitating clump depends upon the value of the 
ratio of the turbulent damping time to the free-fall time 
idamp/ift- If idamp/iff > 1, the turbuleuce fragments the gas 
into self-gravitating cores with an IMF-like mass spectrum. 
These simulations have a star formation efficiency of ~ 45 
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per cent, that is likely an upper limit to the true star forma- 
tion efficiency of the gas. These simulations represent a clus- 
tered mode of star formation. Conversely, if idamp/iff < 1, 
a few large self-gravitating objects arise. The star forma- 
tion efficiency is quite high, ~ 70 per cent. While feedback 
from star formation could possibly reduce this efHciency, it 
is difficult to see how outflows could unbind dense critical 
cores. These simulations are thus more likely to represent 
the conditions of isolated star formation. 

The oblique shocks out of which the cores form impart 
angular momentum to the nascent cores. The magnitude of 
this specific angular momentum is ~ CsL/4, the product 
of the sound speed and the wavelength of the most ener- 
getic turbulent mode. For our simulations, this works out 
to ~ lO^^cm^s"^ - the sam e order of magnitude as seen in 
observations of cloud cores [Goodma n et al.lll993r . but ~ 3 
orders of magnitu de above the values f or the specific angular 
momenta of stars iStassun et al.ll999r . The specific angular 
momentum distributions that we measure for our simulated 
cores peaks roughly at this value of CsL/4, and spans an or- 
der of magnitude in specific angular momentum above and 
below this value. 

T he mass distribution of bou nd cores is well-matched 
by the lPadoan fc Nordlundl i2002l) mass spectrum for cores 
forming out of shocks, in terms of both the width of the 
distribution (related to the current root-mean-square Mach 
number of the fluid motions) and high-mass slope a = —1.29 
(related only to the power spectrum of the turbulence). 
The latter in particular also agrees well with the observed 
Salpeter mass distribution a — —1.35 for molecular cloud 
cores. When scaled to appropriate values, the peak is lo- 
cated at ~ O.lm©, approximately where observations sug- 
gest a peak in the mass spectrum occurs. 

We obtain density profiles for the most massive and dy- 
namically evolved cores in our simulations of p oc r~^. The 
less-developed cores have considerably shallower profiles, as 
significant turbulent motions remain that allow the core to 
support more mass than isothermal pressure along can pro- 
vide. We observe that the collapse proceeds in an aspheri- 
cal manner, preferentially occurring in one direction. This 
is to be expected for a core that is rotating, as rotational 
support will inhibit collapse motions. In a real system, ro- 
tational support can be bypassed through viscous and mag- 
netic forces, but we do not include these effects in our sim- 
ulations. The specific angular momenta of these cores has a 
similar directional departure from spherical symmetry. 

The success of purely hydrodynamic models in deter- 
mining many of the key features of star-forming regions 
brings into question the role of magnetic fields. As mag- 
netic energies tend to be of the same or der of magnitude as 
turbulent energies in molecula r clouds JMvers fc Goodmanl 
Il988l : iBertoldi fc McKeelll992ll . magnetic pressure support 
will likely reduce the SFE of these clumps by making cores 
sub-critical that would otherwise be slightly bound. Mag- 
netic fields will also allow the cores to transport angular 
momentum more efficiently, reducing the total specific an- 
gular momentum of the protostar. Related to this is the fact 
that magnetic fields, coupled with the rotation of the pro- 
toplanetary disk, will drive jets that can inject energy back 
into the cloud. Finally, magnetic fields may serve to increase 
the overall lifetime of the molecular cloud. 
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APPENDIX A: COMPARISON OF 
WATERSHED FLUCTUATION-FINDING 
ALGORITHM WITH CLUMPFIND 

The algorithm we use to find fluctuations, described in Sec- 
tion |5|^ flnds all of the local maxima in the density and all 
of the grid cells whose local gradient vector points towards 
that maxima. Thus, it identifles both gravitationally bound 
cores as well as unbound densit y fluctuations. 

T he CLUMPFIND algorithm JWiUiams. de Geus fc Blitz| 
Il994h is commonly used to find clumps or cores in observa- 
tional maps of star formation regions. This algorithm works 
by finding all of the grid cells with intensity values between 
two contour levels. The algorithm then determines if these 
regions are isolated from any other region or already-found 
core; if the region is isolated, it is labelled a new core. If 
not, the algorithm determines which of the existing cores 
each cell is assigned to. 

We compare the results of these two core-finding codes 
in order to provide a reference point for our watershed algo- 
rithm, to highlight its strengths and weaknesses. We modi- 
fied the CLUMPFIND algorithm to account for a periodic grid, 
by allowing clumps to wrap around the boundaries. We ran 
the CLUMPFIND algorithm on Run B5 (binned to a 128^ grid) 
with a resolution of 0.05 and 0.1 dex in density (so that each 
contour level is 10"°^ = 1.12 or 10°'^ = 1.26 times the previ- 
ous level). For the 0.1 dex resolution search, we find a total 
of 888 fluctuations, of which 259 are bound. These bound 
cores contain a total of 45 per cent of the mass. For the 0.05 
dex resolution search, we find a total of 1319 fluctuations, 
of which 317 are bound (containing 41 per cent of the total 
mass). For comparison, our watershed algorithm found 2886 
fluctuations, of which 207 where bound; 42 per cent of the 
mass is contained in bound cores. The watershed algorithm 
finds significantly more objects that are not destined to col- 
lapse, as it is far more sensitive to small density fiuctuations 
(although only 513 objects are found that have density con- 
trasts between central density and mean surface density that 
is less than 1.25, and only 63 have density contrasts less than 
1.12; none of these are bound. 

The mass distributions for the two clumpfind searches 
are compared with the results of the watershed algorithm 
in Fig. lAll For both clumpfind resolutions, we find good 
agreement with the watershed results for m/mj > 0.002. 
We also find excellent agreement between the distributions 
for the bound cores, which in turn agree very well with 
the Padoan-Nordlund mass distribution. However, the wa- 
tershed algorithm finds significantly more unbound fluctua- 
tions. As our goals are to compare the nature of the bound 
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cores with the unbound fluctuations, we favour the water- 
shed algorithm. 

We have also found that our watershed algorithm exe- 
cutes significantly faster than clumpfind. The clumpfind 
algorithm at 0.1 dex took a comparable amount of time to 
execute as the original simulations did using zeusmp; at 0.05 
dex it took nearly 5 times as long to execute as the zeusmp 
code. The watershed algorithm, for comparison, only took 
about 0.1 times as long. Coupled with the fact that we pre- 
fer the watershed algorithm's definition of a fluctuation, we 
feel that the watershed algorithm is the superior method for 
finding fluctuations and cores in our turbulent simulations. 
However, clumpfind is probably better suited for finding 
cores in observations, as it can more naturally account for 
noise, simply by specifying an appropriate contour level. 
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